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Abstract 

Using data from the FOCUS (E831) experiment at Fermilab, we present new measurements 
for the Cabibbo-suppressed decay mode 7r~7r"'"7r~7r"'". We measure the branching ratio 

r(D° 7r+7r-7r+7r-)/r(Z)° K-tt+tt-tt+) = 0.0914 ± 0.0018 ± 0.0022. An amplitude analysis 
has been performed, a first for this channel, in order to determine the resonant substructure of this 
decay mode. The dominant component is the decay ai(1260)+7r'", accounting for 60% of the 

decay rate. The second most dominant contribution comes from the decay D° /9(770)°p(770)°, 
with a fraction of 25%. We also study the ai(1260) line shape and resonant substructure. Using 
the helicity formalism for the angular distribution of the decay p(770)^ p{770)^ , we measure 

a longitudinal polarization of Pl = (71 ± 4 it 2)%. 

PACS numbers: 13.25Ft,13.30Eg,13.87Fh 
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I. INTRODUCTION 



Hadronic decays of charm mesons are an important tool for understanding the dynamics of 
the strong interaction in the low energy regime. Hadronic decays of D mesons typically have 
a rich resonant structure and a small nonresonant component. Scalar mesons are abundant 
products of three-body decays which have a pair of identical pions in the final state. Due to 
these features, three-body decays of D mesons have been extensively used to study the tttt 
and Kti systems, with eniphasis on the S-wave component of their amplitudes and on the 
scalar resonances 21, 1^, 0, Is]. With D decays one can continuously cover the vrvr and Kn 
mass spectra from threshold up to 1.7 GeV/c^ {Mr, — rriT^), filling the existing gaps between 
the 6] and the CERN-Miinich data [3], in the tttt case, and between threshold and 825 
MeV/c^, where the LASS data g] on Kit scattering start. Another interesting feature of 
D decays is that the bulk of the hadronic decay width can be described in terms of simple 
tree-level quark diagrams. There seems to be, at least for the vector, axial-vector and tensor 
mesons, a strong correlation between the final state quarks — the spectator valence quark 
plus the ones resulting from the weak decay of the c quark — and the quark content of the 
observed resonances in the intermediate states. This connection allows some insights on the 
nature of light mesons. The information provided by D decays is, therefore, complementary 
to that of traditional scattering experiments. 

Much less information is available concerning the tttttt / Kmr systems and the axial- vector 
resonances. The picture emerging from many different Dalitz plot analyses of D decays 
reveals a well defined pattern: final states that can be associated with a simple spectator 
amplitude (PF-radiation), in which the virtual W couples to a vector or axial- vector meson, 
have large branching fractions compared to final states in which the W is coupled to a 
pseudo-scalar meson. Following this pattern, one expects the ai(1260)+7r^ channel to be 
the dominant resonant component of the 7r^7r+7r^7r+ decay (charge conjugation is 

always implied, unless stated otherwise). This decay, therefore, can be used to study the 
resonant substructure of the axial- vector meson ai(1260), as well as its line shape, from 
threshold up to 1.72 GeV/c^. 

In the study of light mesons from D decays there is one potential difficulty: in hadronic 
decays the tttt and Kn pairs, as well as the -k-r-k / Kti-k systems, are part of a few-body 
strongly interacting system. In principle one must account for final state interactions (FSI) 
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between all decay particles. It is still an open question whether the FSI are strong enough in 
three-body decays to distort the pure tttt/J^tt scattering amplitudes. In any case, the effects 
of the FSI should become more important as the number of final state particles increases. 

Most amplitude analyses employ the so called isobar model: a coherent sum of resonant 
amplitudes weighted by constant complex coefficients. The constant phases account for re- 
scattering effects between the isobar and the other decay particles. This approach to FSI 
may be sufficient in some cases, but may be too simplistic in four-body decays. In the most 
general case, the FSI would depend on energy, with a smooth variation across the phase 
space. A correction to the isobar model would then be necessary in order to incorporate 
the energy dependent FSI. As we will show, the energy dependent FSI can probably not be 
ignored in four-body decays of D mesons. 

There is an additional motivation for the study of four-body D meson decays. Charmless 



decays of B mesons are a promising tool for the study of CP violation. 
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n particular, the 



12|. The B 



decays B pp have been used to extract the CKM angle a [ol, |lo|, 
ai(1260)7r mode, however, also leads to the same innni final state. The B^ —>■ ai(1260)^7r" 
channel has a branching fraction which is one order of magnitude larger than that of B'^ — 



p'p' 



13 



14| . There is no measurement of the B^ — > 01(1260)^71*^ branching ratio, but since 



this decay proceeds mainly via the same tree-level W radiation amplitude as that of the 
p^p*^ channel, one expects these two modes to have comparable rates. The B —>■ ai(1260)7r 
decay, thus, accounts for a large fraction of the B titititc decay, and a full amplitude 
analysis would be necessary in order to isolate the B pp contribution, as in the case of 
the —>■ 71^71^71^71^ decay. All the systematics — form factors, resonance line shapes, 
representation of the S-wave components, angular distributions, FSI, etc. — are common to 
both decays. One can think of the 7r~7r+7r~7r+ decay as a prototype of the B — > TrTrTrvr 

decay. 

In this paper, we present a new measurement of the relative branching ratio (BR) r(Z)° — >■ 
7r^7r+7r^7r+)/r(Z}*^ —>■ K^7t^7i^7i^) using data from the FOCUS experiment. For the first 
time an amplitude analysis has been performed to determine the 7r^7r+7r^7r+ resonant 

substructure. 
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II. THE FOCUS EXPERIMENT 



FOCUS, an upgraded version of E687 [l5|, is a charm photo-production experiment which 
collected data during the 1996-97 fixed target run at Fermilab. Electron and positron beams 
(typically with 300 GeV endpoint energy) obtained from the 800 GeV Tevatron proton beam 
produce, by means of bremsstrahlung, a photon beam which interacts with a segmented BeO 
target [16[. The mean photon energy for reconstructed charm events is ~ 180 GeV. A system 
of three multi-cell threshold Cerenkov counters performs the charged particle identification, 
separating kaons from pions up to a momentum of 60 GeV/c. Two systems of silicon 
micro-vertex detectors are used to track particles: the first system consists of 4 planes of 
micro-strips interleaved with the experimental target [17t] and the second system consists of 
12 planes of micro-strips located downstream of the target. These high resolution detectors 
allow the identification and separation of charm primary (production) and secondary (decay) 
vertices. The charged particle momentum is determined by measuring the defiections in two 
magnets of opposite polarity through five stations of multi-wire proportional chambers. 



III. THE D° TT 7r+7r-7r+ SAMPLE 

The final states are selected using a candidate driven vertex algorithm [l^ . To minimize 
systematic errors on the measurements of the branching ratio, we use identical vertex cuts 
on the signal and normalizing mode. 

Secondary vertices are formed from the four candidate tracks. The momentum of the 
resultant candidate is used as a seed track to intersect other reconstructed tracks and 
to search for a primary vertex. The confidence levels of both vertices are required to be 
greater than 1%. The mean energy of the candidates in the LAB frame is 85GeV. This 
is equivalent to an average decay length of 1cm, allowing a good spatial separation between 
the production and decay vertices. Once the position of the production and decay of the 

candidate is determined, the distance L between the vertices and its uncertainty 
are computed. The ratio L/ai is the most important variable for separating charm events 
from non-charm prompt backgrounds. Signal quality is further enhanced by cutting on 
the isolation variables, Isol and Iso2. The isolation variable Isol requires that the tracks 
forming the D candidate vertex have a confidence level smaller than the cut to form a vertex 
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with the tracks from the primary vertex. The Iso2 variable requires that all remaining tracks 
not assigned to the primary and secondary vertex have a confidence level smaller than the 
cut to form a vertex with the D candidate daughters. In addition, we require the secondary 
vertex to lie outside of the segmented targets (Out- of- Material cut, or OoM), in order to 
reduce contamination due to secondary interactions. The OoM variable is actually the the 
distance between the secondary vertex and the edge of the nearest target segment divided by 
the uncertainty on the secondary vertex location. We have also applied a cut on vtXgcore, a 
variable built from the OoM, L / and the confidence level of the secondary vertex, in order 
to explore the correlations between these variables. Each of these variables is normalized to 
its maximum value, so vtXscore ranges from to 1. This cut allows us to further reduce the 
background without applying tighter cuts on each of these variables. The set of cuts that 
yield the h^ti^ti^ti^ signal with the best statistical significance is L / ai > 10, Isol 

and Iso2 < 10%, vtXscore > 0.15 and OoM > 1. 

The only difference in the selection criteria between the Ti^n^n^n^ and — >■ 

K^7r+7r^7r+ decay modes lies in the particle identification cuts. The Cerenkov identification 
cuts used in FOCUS are based on likelihood ratios between the various particle identification 
hypotheses. These likelihoods are computed for a given track from the observed firing 
response (on or off) of all the cells that are within the track's {f3 = 1) Cerenkov cone for 
each of our three Cerenkov counters. The product of all firing probabilities for all the cells 
within the three Cerenkov cones produces a x^-like variable Wj = —2 log(Likelihood) where i 



ranges over the electron, pion, kaon, and proton hypotheses ISj . The kaon track is required 



to have Ak = W-,, — Wk > 3, and all pion tracks are required to be separated by less than 
5 from the best hypothesis, that is = Wmm — W^tt > —5. 

Using the set of selection cuts just described, we obtain the invariant mass distribution for 
7T~ TT~^ 7T~ 71^ showu in Fig. [H Although the Cerenkov cuts considerably reduce the reflection 
peak (from — ^ K^ti^ti^ti^) to the left of the signal peak, there is still a large distortion 
of the background due to this surviving contamination. The tt^h^ti^ti^ mass plot is fit 
with a function that includes two Gaussians with the same mean but different sigmas to 
take into account the variation in momentum resolution of our spectrometer as a function 
of particle momentum 15|, and two exponential functions for the combinatorial background 
and for the K~tt~^7t^7i^ reflection. A log-likelihood flt gives a signal of 6360 ±115 

7r^7r+7r^7r+ events, over a background of 769 ± 15 events. 
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FIG. 1: Invariant mass distributions used to determine the ratio of branching fractions. The upper 
plot is the 7r~7r"'"7r~7r+ signal; the lower plot is the normalizing channel i^~7r+7r~7r+. The fits (solid 
curves) are explained in the text; the dashed line shows the background. The hatched area on the 
7r~7r+7r~7r+ signal corresponds to the events used in the amplitude analysis. 

The large statistics K^n'n'^n'^ mass plot is also fit with two Gaussians with the same 
mean and different sigmas plus a second-order polynomial for the background. A log- 
hkelihood fit gives a signal of 54156 ± 267 D° — > K~7r~7r'^7r'^ events. 



IV. RELATIVE BRANCHING RATIO 



The evaluation of relative branching ratios requires yields from the fits to be corrected for 
detector acceptance and efficiency. These differ among the various decay modes, depending 
on the Q values and Cerenkov identification efficiency. 

Prom the Monte Carlo simulations, we compute the relative efficiencies (with statistical 
error only): 

e(D° X~7r+7r""7r+) 

\ J = 0.7891 ± 0.0004. (1) 

Using the previous results, we obtain the following value for the branching ratio: 



8 



Systematic uncertainties on the branching ratio measurement may come from different 
sources. We estimate the systematic uncertainty on the yields and on the efficiencies consid- 
ering three independent types of contributions: the split sample component, the fit variant 
component, and the component due to the particular choice of the vertex and Cerenkov 
cuts. 

The split sample component takes into account the systematics introduced by a residual 
difference between data and Monte Carlo, due to a possible mismatch in the reproduction 
of the D'^ momentum and the changing experimental conditions of the spectrometer during 
data collection. This component has been determined by splitting data into four independent 
subsamples, according to the momentum range (high and low momentum) and the 
configuration of the vertex detector, that is, before and after the insertion of the upstream 
silicon system. A technique, employed in FOCUS and modeled after the S-factor method 
from the Particle Data Group [19], was used to try to separate true systematic variations 
from statistical fluctuations. The branching ratio is evaluated for each of the four statistically 
independent subsamples and a scaled error a (that is the errors are boosted when / {N — 
1) > 1) is calculated. The split sample error agput is defined as the difference in quadrature 
between the statistical error returned by the fit on the unsplit data set, agtat, and the scaled 
error, if the scaled error exceeds the statistical error, asput = — oltat- 

Another possible source of systematic uncertainty is the fit variant. This component 
is computed by varying, in a reasonable manner, the fitting conditions for the whole data 
set. In our study we fixed the widths of the Gaussian to the values obtained by the Monte 
Carlo simulation, we changed the background parametrization (varying the degree of the 
polynomial), and we use one Gaussian instead of two. The latter is the dominant contri- 
bution to the fit variant error. Finally the variation of the computed efficiencies, both for 

—>■ 7r^7r+7r^7r+ and the normalizing decay mode, due to the different resonant substruc- 
ture simulated in the Monte Carlo has been taken into account. The BR values obtained by 
these variants are all a priori equally likely, therefore this uncertainty can be estimated by 
the r.m.s. of the measurements. 

We estimate the cut component systematic error by varying vertex and particle identifi- 
cation cuts one at a time. We varied the confidence level of the secondary vertex from 1% 



TABLE I: Contributions to the systematic uncertainties of the branching ratio 



Source 


Systematic error 




Spht sample 
Fit Variant 
Set of cuts 


0.0010 
0.0012 
0.0016 




Total systematic error 


0.0022 




TABLE IL Comparison with other experiments. 


Experiment 


r(Z)0 ^ 7r-7r+7r-7r+)/r(Z)° ^ K-tt+tt-tt+) 


Events 


FOCUS (this result) 
CLEO-c [20] 
BES [2]J 
E687 [22] 


0.0914 ±0.0018 ±0.0022 
0.097 ±0.002 ±0.003 
0.079 ± 0.018 ± 0.005 
0.095 ±0.007 ±0.002 


6360 ± 115 
7331 ± 130 
162 ± 20 
814 ± 26 



to 50%, Isol and Iso2 from 10~^ to 1, L / ai from 6 to 20, from —6 to —2, and vtxgcore 
from to 0.2. Analogously to the fit variant, the cut component is estimated using the 
standard deviation of the several sets of cuts. Actually, this is an overestimate of the cut 
component because the statistics of the cut samples are different. 

Finally, adding in quadrature the three components, we get the final systematic errors 
which are summarized in Table [H 

The final result is shown in Table HTl along with a comparison with previous measurements. 

V. AMPLITUDE ANALYSIS 

A fully coherent amplitude analysis was performed in order to determine the resonant 
substructure of the 71^71+71^71+ decay. This is, to our knowledge, the first such 

measurement for this channel. The amplitude analysis was performed in the framework of 
the isobar model on a sample of 6153 events, corresponding to the hatched area in Fig. [1] 
(events with m^.,^ within 20 MeV/c^ of the mass). 
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A. A model for the —>■ vr tt+vt Tr"*" decay 

The decay mode ti~ti^ti~ti^ is Cabibbo-suppressed and may proceed through many 

intermediate resonant states. Considering only tree-level amplitudes, possible contributions 
come from two-body decays such as D ^ ai(1260)+7r-, D p{770f p{770f , D p{770fS 
and D SS, where S* is a scalar meson {S = /o(600) or cr, /o(980), /o(1370)). The 
7r^7r+7r^7r+ can also result from three-body nonresonant decays D Rtttt {R = /o(600) or 
a, p(770)°, /o(980), /2(1270), /o(1370), p(1450)°). Including all possible contributions, one 
could have over twenty different resonant amplitudes leading to the n^n^Ti^ii^ final state. 

In addition to the resonant modes, the 7r^7r+7r^7r+ final state could also result from a 
four-body nonresonant decay. Usually the nonresonant component is assumed to be uniform, 
)ut this may not be a reasonable assumption even in the simpler case of three-body decays 



23l |. To our knowledge there is no phenomenological model for a non-uniform nonresonant 
amplitude in four-body decays. 

The existence of many possible intermediate states leading to the 7t~tt~^tt~7c~^ final state 
makes the amplitude analysis of this decay very challenging. Many of these intermediate 
states involve broad resonances that populate the whole phase space. A model having a large 
number of overlapping amplitudes gives rise to large interference terms which are difficult to 
control in a five- dimensional space. The presence of these interference terms is easily detected 
when the sum of decay fractions greatly exceeds 100%. In such cases many local minima 
exist in the parameter space, with very similar likelihood values. With such a model one 
may find a mathematical solution to the fit problem, losing, however, the physical meaning 
of it. An additional difficulty comes from the fact that in the 7r^7r+7r^7r+ final state there 
are two pairs of identical pions and therefore amplitudes must be Bose-symmetrized. In the 
chain ai(1260)+7r^, ai(1260)+ p'^^r^, tt+tt^, for instance, one can combine 

the four pions in eight possible ways. Angular distributions, a clean signature of specific 
modes, tend to be smeared out by the Bose-symmetrization. This could be minimized using 
flavor tagging, since a large fraction of the comes from the decay D*^ D^{D^)7i^. 
In our case that would reduce the sample size by a factor of five. Finally, as the number 
of final state particles increases the more important the final state interactions tend to be. 
The FSI may play a significant role in the —>■ 7r^7r+7r^7r+ decay. Unfortunately this is a 
very difficult problem to deal with, even with a phenomenological approach. 
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The strategy adopted in this analysis is to start with a set of amphtudes corresponding 
to modes that are expected to be dominant. Once a stable solution is achieved, modes with 
marginal contributions are replaced by other ones, until a final set is reached. 

We built a baseline mode guided by the n^n^ projections, by MC simulation of the 
different channels and by our own experience. The 7r^7r+ mass projections are shown in Fig. 
[21 The plot in Fig. [2]^a) has four entries per event, corresponding to the four different vr^vr"'" 
combinations one could form in each event. We adopt the particle labeling vrf 7r^7r3"7r^, 

so the identical particles are ni/n^ and 112/ Hi- In our sample we do not distinguish between 

and D^, and the like charge pions are randomized to avoid any ordering with respect to 
their momenta. Consequently, all four pions have the same momentum distribution, making 
the four 7r^7r+ mass projections indistinguishable. 

One can form two combinations of 7r^7r+ pairs, vrf vr^/TTg^Tr^'" or Trjf Tr^'^/Trg^vr^. One can 
plot the highest and lowest tt'tt'^ invariant mass for each pair combination. These are shown 
in Fig. [2t^b)-(c), respectively. In these plots there are two entries per event. 

The p(770)° is clearly the dominant tt^tt^ resonance, although other resonances, like 
the a, /o(980) and /2(1270) may also be present, even without leaving a clear signature in 
the vr^TT^ mass projection. The p(770)'' signal could originate from the decay of the axial- 
vector meson ai(1260), from the decay of the type P^P^, P^Ri or even from 
£)0 _^ pO^+^-^ Figure 2(d) shows the tt^tt^ projection of Monte Carlo (MC) simulations of 
the decay ai(1260)+7r^, with the Oi decaying to p°7r+ (S- and D-wave) and to an'^. 

Our baseline model, therefore, includes contributions of three types. The first type is 
the —>■ ai(1260)^7r~ chain. For the ai(1260) resonant substructure, we consider three 
channels: ai — > p(770)'^7r^, the p^n^ being in a dominant S-wave state with a small D-wave 
component, plus cttt^ . The second type of contribution is the —>■ p^p^ decay, in 

three possible helicity states. Finally, we consider decays of the type Rtt^tt~, with 

i? = a,pO,/o(980) and/2(1270). 

B. Formalism 

The formalism used in this amplitude analysis is a straightforward extension to four- 
body decays of the usual Dalitz plot fit technique. The —>■ 7r^7r+7r^7r+ signal amplitude 
is represented by a coherent sum of individual amplitudes, each corresponding to a possible 
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FIG. 2: Invariant mass distribution for vr^vr"'". (a) All four 7r~7r"'' combinations; (b) the TT~7r~^ 
combinations with highest mass; (c) the 7r~7r+ combinations with lowest mass; (d) MC simulation 
of the 7r~7r+ mass distribution for D° — ai(1260)''"7r'~: the solid line corresponds to the 
/9(770)^7r+ S-wave, the dashed line is the af — p(770)°7r+ D-wave, and the dotted line is — (T7r+. 

intermediate state and weighted by constant complex coefficients Ck, A — '^CkA/.. The 
background amplitude is represented by an incoherent sum of amplitudes corresponding to 
the different types of background events. 

A likelihood function is built from the signal and background amplitudes, incorporating 
detector resolution effects and acceptance correction. The optimum set of constants c^, 
the only unknowns, is obtained by an unbinned maximum likelihood fit. The background 
amplitude is kept fixed during the fit procedure. 



1. Phase space 

The kinematics of the 7r~7r+7r~7r+ decay is unambiguously determined by five 

degrees of freedom. Wc choose a set of five two-body invariant masses squared — the 
four different TT^Tr"*" combinations, S12, su, S23, S34 plus S13 — to define an event in the five- 
dimensional phase space. Any kinematic variable, like acoplanarity and helicity angles or 
three-body masses, can be easily expressed in terms of this set of invariants. Unlike the case 
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of three-body decays, the phase space density is no longer constant. 



2. Signal amplitudes 

The amphtudes for the intermediate states, Aj, are phenomenological objects. Each signal 
amplitude is built as a product of relativistic Breit-Wigner functions (BW), an overall spin 
amplitude, M., accounting for angular momentum conservation at each decay vertex, and 
form factors, Fi {I is the orbital angular momentum of the decay vertex), accounting for 
the finite hadron size. This is the standard structure for resonant amplitudes (see note 



of Dalitz plot formalism in 19|]). The amplitude for the decay chain ai(1260)^7r , 

ai(1260)+ ^ p°7r+ (S-wave), p° ^ tt+tt-, is 

A = F^ ^ BW"i X X BW^ x Ff x A^. (3) 
The amplitude for the decay chain P1P25 Pi ~^ tt+tt^, p2 ~^ tt+tt^ is 

A = Ff ^ BW^^ X X BW^2 x Ff x A^, (4) 
and for the chain /2(1270)7r+7r^, /2(1270) tt+tt^ the amplitude is 

^ = X BW^^ X fI"" xM. (5) 

All signal amplitudes are Bose-symmetrized due to the existence of two pairs of identical 
particles in the final state. 

n 

We use the Blatt-Weisskopf damping factors [24] as form factors for vertices involving 
spin-1 and spin-2 resonances. These form factors have one free parameter, r/j, which is 
related to the resonance size. All form factors are slowly varying functions of energy. As a 
systematic check we varied the value of r^j. The fit result is not very sensitive to the value 
of this parameter in the range 1 < r^j < 5 GeV~^. The form factor parameter is fixed at 
tr = 3.0 GeV"^ 

For the spin amplitudes we use the Lorentz invariant amplitudes for the decay sequences 
afn^ and D — > Rtt^tt^ {R being a vector or tensor resonance). These Lorentz 



invariant amplitudes are built using the relativistic tensor formalism 25|, 1261] which combines 



particle 4-momenta with polarization vectors in a rotationally invariant way. In vertices 
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involving only the strong interactions, parity conservation imposes an additional constraint. 
For the decay p^p^ the spin amplitudes are written using the helicity formalism. 

Explicit formulae for the angular distributions M are presented in the Appendix. 

The relativistic Breit-Wigner representing the /2(1270) and a resonances has an energy 
dependent width, with the approximation r^^(s) = rtot(s), 

s- So + iy/s^T totis) ' 

where 

In the above equations s is the vr+vr" mass squared, p* is the break-up momentum of the 
resonance, Sq the resonance nominal mass squared, and = p*{sq). 

The p° line shape is better described when the interference with the C(j(782) is included. In 
spite of the tiny branching fraction of the ti;(782) —>■ tt~^7t~ , the effect of the p-u interference 
is remarkable and clearly distorts the pure p° line shape, as one can see by comparing the 



TT+TT with the TT^Tc'^ mass distributions from 



given by the Crystal Barrel Collaboration 29 1 



In this analysis we use the line shape 



For the /o(980) we used the Flatte formula [30] of a coupled channel Breit-Wigner 
function. The /o(980) parameters — gn = 0.20 ± 0.04, and qk = 0.50 ± 0.20 and 
mo = (0.957 ± 0.008) GeV/c^ — are obtained by a fit to the FOCUS D+ 71+71-71+ 
Dalitz plot, where the /o(980)7r"*" is the dominant component. 

The /2(1270) Breit-Wigner parameters are fixed at the PDG values 13], whereas for the 
a we used the E791 parameters It is well known that the a pole cannot be obtained 

n 

from a simple Breit-Wigner [27|. Nevertheless, good fits have been obtained by different 
experiments, in different channels, using Breit-Wigners with similar parameters. There is 
no consensus on the correct way to parameterize this low mass, broad scalar state. The focus 
here, however, is not on the a properties, but on the ai(1260) line shape, which depends on 
its resonant substructure. For this purpose it suffices to use the Breit-Wigner formula as 
an effective representation of a scalar component of the ai(1260) resonant substructure. In 
Section VI. D we discuss the sensitivity of the fit results to the values of the a parameters. 
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3. The ai(1260) line shape 



One of the main goals of this analysis is to determine the ai(1260) line shape and its 
resonant substructure. These are correlated: we need to know the resonant substructure 
in order to determine the line shape because the resonant substructure defines how the 
ai(1260) width, r^ot(s) depends on the ti^ti^ti^ mass squared s. Given the functional form 
of rtot(s), we represent the ai line shape by the same Breit-Wigner formula of Eq. ([6]). The 
total width of the ai(1260) is given by 

We assume that r2^07r+('^) ~ ^tl+n-^^^ '^^^ '^i '^^^ decay to KKn via /o(980)7r 
or K*{892)K. The opening of these channels introduces a cusp-like effect in the ai width. 
Since no evidence of the mode /o(980)7r has been reported so far, we assume the KKir 
partial width is entirely due to K*{892)K. The value of the coupling constant gx^x taken 
from the analysis of the decay r tt'^tt^tt+z/^ by CLEO [31]. The coupling constant Qx^x 
was also measured by the Belle Collaboration [s^, which found a value three to five times 
larger than that measured by CLEO. The uncertainty on the value of g^^x ^ source of 
systematic uncertainty. 

The partial width r'^'^+^_{s) is obtained by integrating the tt+tt"??"'" Dalitz plot, 

TZ+^S) « ds^dS2\AaM,S2)\\ (9) 

with si, S2 being the invariant mass squared of the two possible vr+vr^ combinations and 

3 

S2) = ^CiA(si, S2), (10) 

i=l 

where q are complex coefficients and Ai{si,S2) the amplitudes for each of the ai decay 
modes. We fix the partial width by demanding that ^2l-+^~{s = Sq^) = Fq^ with Sq^ and 
Fq^ being real input parameters to be determined. 

In practice, we used an iterative procedure: assuming initial values for Sq^ and Fq^, and 
for the constants q, the total width Ftot(s) is computed. Having computed Ftot(s), we 
perform the five-dimensional fit, which returns new values for the constants q; the function 
^tlr+TT' ('^) then recomputed using the same values of Sq^ and Fq^ and the new values of the 
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constants Cj. The process converges after a few iterations. The whole procedure is repeated 
many times, scanning over the values of Sa^ and Fq^ until we find the combination that 
optimizes the fit result. 

4- Background 

The background composition was determined by inspection of the data on the side band 
to the right of the signal in the 7r~7r"'"7r~7r"'" mass spectrum. We consider two types 
of background events: random combinations of a p(770)° and a tt^tt"*" pair and random 
combinations of 7r~7r+7r~7r+. We estimate a fraction of 80% for random combinations of 
7r~7r+7r~7r+ and 20% for p(770)° plus a 7r~7r+ pair. Random combinations of two p° might 
be also present, although with a small contribution. The background fractions, as well as 
its composition, is a source of systematic uncertainty. 

We assume the random 7r~7r"'"7r~7r"'" combinations to be uniformly distributed in phase 
space, whereas for the other background we assume a Breit-Wigner with no form factors and 
no angular distribution. The overall background distribution is a weighted, incoherent sum 
of the two components described above. The relative background fractions, bk, are fixed in 
the fit. 

5. Detector resolution and acceptance correction 

The finite momentum resolution causes a smearing of the phase space boundary. It is 
easier to understand this effect if we consider the case of a three-body decay and the Dalitz 
plot. In a sample where the three-body masses have only one well defined value, the Dalitz 
plot boundary is also well defined and unique. But in a sample where we have a distribution 
of three-body masses, G{M), the Dalitz plot boundary is no longer well defined but, rather, 
it is a superposition of boundaries, each one defined by its corresponding value of the three- 
body mass. Dividing the three-body mass distribution into fine bins, one can think of the 
observed Si2 x S23 distribution as a superposition of Dalitz plots, each one weighted by 
G{Mi). The same reasoning applies to four-body decays. We incorporate the effect of the 
momentum resolution, both in the signal and background probability distribution functions 
(PDF), by weighting each event by the Att mass distribution. With this approach, the signal 
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and background PDF depend on the five invariants and on the An mass as well. 

The acceptance correction is also applied on an event-by-event basis. The acceptance is a 
function of the five invariants. In the case of the amplitude analysis, the absolute value of the 
acceptance function is irrelevant, since it implies an overall constant factor multiplying the 
likelihood. All one really needs to know is the acceptance at a given phase space cell relative 
to its neighbor cells. We use a five-dimensional matrix for the acceptance correction rather 
than a continuous function. We used the full FOCUS Monte Carlo simulation (ROGUE) to 
generate a very large sample of D'^ — > 7r~7r+7r~7r''" events with a constant matrix element. 
We applied the same selection criteria to this sample as those used for real data. A mini-MC 
simulation of the 7r~7r+7r~7r+ was used to generate the phase space distribution. Each 

of the five axes were divided into bins. Two five-dimensional arrays were filled, one with the 
events of the full simulation that passed all cuts and the other with the phase space events. 
The acceptance matrix was formed by dividing the number of ROGUE events in each cell 
by the number of phase space events in that cell. The acceptance is nearly uniform across 
the phase space, with a 10-15% decrease close to the edges. 

6. Normalization 

The maximum likelihood fit technique requires the likelihood function to be normalized. 
This means one normalization constant for the signal and another for the background prob- 
ability distribution functions. The normalization constant for the background distribution 
needs to be calculated just once, since the background amplitudes and the relative frac- 
tions are fixed during the fit. The normalization constant for the signal, however, must be 
computed at each step of the minimization procedure, since it involves the fit parameters. 

Although it is not necessary, the amplitudes for each intermediate state are normalized 
to unity. This is intended to give the magnitudes Ck a direct physical meaning: the decay 
fraction of each mode is directly proportional to |cfep. 

The normalization is a crucial step in amplitude analysis. The evaluation of the normal- 
ization constants requires a number of phase space integrals. Fortunately this integration 
can be performed just once, since all input parameters in the amplitudes — mass and width 
of resonances, form factors parameter — are fixed during the minimization process. 

We define G(m47r) and 6(m4,r) as functions representing the Att mass signal and back- 
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ground distributions, respectively. The function G(m47r) is a sum of two Gaussians, whereas 
6(m47r) is a sum of two exponentials. Both G(m4^) and &(m47r) were obtained by the fit 
to the Att mass spectrum shown in Fig. [TJ If we define as a point in phase space 
{d(f) = (isi2(isi4(is23(is34(isi3), p(0) as the phase space density, and 6{(j)) as the efficiency, 
and Ak{(j)), Bj{(j)) the unormalized amplitudes, the normalization integrals are 

iVf = J dm^^ G(m4.) J d(t> p{<j)) \Ak{<j))\^ (11) 



and 



Nf = j dm^^ b{mi^) J dcp p{(p)Bj{(p). (12) 
The overall signal and background normalization integrals are 

Ns = j dm^^ G{mi^) j c/^ p(0)e(0) |^ CfcAfe(0) (13) 

and 

Nb = j dm^^ b{mi^) j d^ p{(j))e{(l))^bjBj{(f)), (14) 

with Afc(0) = Afc(0)/iVf and 5,(0) = Bj{(l))/Nf. 

We performed the phase space integration using the Monte Carlo. The m4,r interval is 
divided into bins. For each m^T^ bin a very large sample of MC events, generated with a 
constant matrix element, is used to compute average value of the integrands. These average 
values are then multiplied by the phase space volume and weighted by (7(7714^) /6(m4^). 



7. The likelihood function 

The overall signal and background amplitudes are corrected on an event-by-event basis for 
the acceptance, which is nearly constant across the phase space, and for the finite detector 
resolution, taken into account by multiplying the overall signal distribution by a Gaussian 
factor, G'(m4^), and the background distribution by an exponential function, 6(7714^). The 
normalized signal probability distribution is, thus. 



SpBF{(p]m4^) = j;r^i^4^)u2'^kAk{(j)) , (15) 



Ns 
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and 



Ps((/);m4^) = £(0)p(0)5'pDF(0;"^47r)- 

The normalized background probability distribution is 



(16) 



Si 



PDF(0;m4,r) = -^b{m4^)^bkBk{(l)), 

iVR 



(17) 



and 



Pb(0; m4,r) = £(0)p(0)-BpDf(0; "7,47r)- 



(18) 



An unbinned maximum likelihood fit was performed, minimizing the quantity w = 
— 21og(>C). The likelihood function, C, is 

£= n [Ps{<l>';ml,) + PB{(l>';ml,)]= J] [£(0>(</'O('5pdf(</''; + i^PDF(0'; m^)] 

events events 

(19) 

When we take the logarithm of the likehhood we have 



iogc=j2 log H<Pip{<P')] + Yl [SpM<p'; <j + i?PDF(0^ m^j] (20) 



events 



events 



In the fit parameter space the term X] log [^(0*)p(0*)] is a constant and does not affect 
the position of the minimum. In practice, thus, we minimize the quantity 



w = -2 



—G{m^^) V CfcAfe(0^) + —b{m4n) Y] 6^5^(0*) 



(21) 



The acceptance correction and the phase space density enter only in the overall normal- 
ization constants Ns and Nb- 

Decay fractions are obtained from the coefficients c^, determined by the fit, and after 
integrating the overall signal amplitude over the phase space at m47r — moo, 



fk 



2 ' 



(22) 
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TABLE III: Results from the best fit. The first error is statistical, and the second one is systematic. 



mode 




magnitude 




phase (degrees) 




fraction I 




a^7r~, ai — > p^ir^ (S-wave) 




1. 


(nxed) 




(fixed) 


43.3 


zb 


2.5 zb 


1.9 


a^vr^, ai p^n^ (D-wave) 


0.241 


lb 


0.033 ± 


0.024 


82 lb 5 lb 4 


2.5 


zb 


0.5 zb 


0.4 


a^TT^, ai avr^ 


0.439 


1 

± 


O.OZD ± 


0.021 


193 zb 4 zb 4 


8.3 


1 

± 


0.7 ± 


0.6 


ati:^ (all) 












60.0 


zb 


3.0 zb 


2.4 


p'^P*' (parallel) 


0.157 




0.027 lb 


0.020 


120 lb 7 lb 8 


1.1 


zb 


0.3 zb 


0.3 


p'^ (perpendicular) 


0.384 




0.020 lb 


0.015 


163 zb 3 zb 3 


6.4 


zb 


0.6 zb 


0.5 


p'^ p^ (longitudinal) 


0.624 




0.023 lb 


0.015 


357 zb 3 zb 3 


16.8 


± 


1.0 zb 


0.8 


p^p^ (all) 












24.5 


zb 


1.3 zb 


1.0 


/o(980)7r+7r- 


0.233 




0.019 lb 


0.015 


261 zb 7 zb 4 


2.4 


zb 


0.5 zb 


0.4 


/2(1270)7r+7r- 


0.338 




0.021 lb 


0.016 


317 zb 4 zb 4 


4.9 


zb 


0.6 zb 


0.5 




0.432 




0.027 lb 


0.022 


254 zb 4 zb 5 


8.2 


zb 


0.9 zb 


0.7 


iivr+vr- (all) 












20.0 


zb 


1.2 zb 


1.0 



since the individual amplitudes ^4^(0) are normalized to unity. Errors on the fractions 
include errors on both magnitudes and phases, and are computed using the full covariance 
matrix. 

VI. RESULTS FROM AMPLITUDE ANALYSIS 

The technique described in Section V is applied to the events in the hatched area of Fig. 
1. The decay Oj'^tt^, p^n^ (S-wave) is taken as the reference mode, fixing the 

phase convention and the relative magnitudes of the other contributions. The results of the 
best fit to the data are summarized in Table IIIII The systematic errors are discussed in 
Section VI. E. 

In Fig. [3] the vrvr mass projections from data (diamonds with error bars) is plotted with 
the fit result overlaid (solid histograms). The fit projections are made from a large MC 
simulated sample, including the signal and background PDF's, defined in Eq. ( fTSl) and (fTTj) . 
as well as the reconstruction efficiency and detector resolution. In Fig. [3]^a) the tt^tt^ mass 
is shown. Events in Fig. [3](b) and (c) are the tt^tt^ combinations with highest and lowest 
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FIG. 3: Comparison between data and the result of the best fit. In (a) all four vr^vr^ combinations 
are added. Entries from plot (a) are split into combinations with highest and lowest vr^vr^ mass: 
in (b) we show the Tr^vr^ combinations with highest mass, whereas in (c) we show the 7r~7r+ 
combinations with lowest mass; finally, in (d) we show the 7r~7r^/7r+7r^ mass distribution. In all 
plots the solid histogram is a projection of the fit. 

mass, respectively. In Fig. [3t^d) we show the -n^-n^ j-n^-n^ invariant mass. 

The TT^TT^TT^ mass spectrum is shown in Fig. HI This plot has four entries per event. 
The diamonds with error bars are the data distribution, whereas the solid histogram is the 
fit projection. Note that this distribution does not reflect the pure ai(1260) line shape, for, 
in addition to the Bose-symmetrization, there are contributions from the other modes. 

A further comparison between the fit result and the data can be made with the pion 
momentum distribution, computed in the D rest frame. Recall that we do not distinguish 
between D° and -D°, and that there is no particle ordering according to its momentum. 
This means that all pions must have the same momentum distribution and that we can add 
them into a single plot. The pion momentum distribution is shown in Fig. ([5]), where the 
diamonds represent the data points, and the solid histogram is a MC simulation. 
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FIG. 4: The vr+vr 7r+/7r tt+tt mass spectrum. A projection of the fit (soUd histogram) is overlaid 
with the data points (diamonds with error bars). 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



pion momentum in D rest frame (GeV) 



FIG. 5: The pion momentum distribution in the D rest frame (four entries per event). The 
momentum distribution for all pions should be the same, so in this plot there are four entries per 
event. A projection of the fit (solid histogram) is overlaid with the data points (diamonds with 
error bars). 
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FIG. 6: The cosine of the acoplanarity angle distribution. The contribution of the p^p^ is 

only 25%, so the original polarization of the two vector states is diluted. 

We can also compute the acoplanarity angle. In the D RR decay this is the angle, 
measured in the D rest frame, between the planes defined by the decay particles of each of 
the resonances. Note that in the 7r^7r+7r^7r+ decay the p^p^ mode accounts for only 

one fourth of the decay rate, so the distribution of the acoplanarity angle no longer carries 
the information on the polarization state of the two vector mesons. The distribution cosine 
of the acoplanarity angle is shown in Fig. [61 Once more, the diamonds represent the data 
points, and the solid histogram is a projection of the fit. 

A. Goodness-of-fit 

The TT+vr" and the n~^7i^7i^ mass spectra, as well as the pion momentum and the cosine of 
the acoplanarity angle distribution, are reasonably well represented by our model. These are, 
however, inclusive distributions, indicating that the model represents the general features of 
the data. Any quantitative statements about the fit quality should be made by inspecting 
the five-dimensional phase space. 

The goodness-of-fit is assessed using a x^-like test. The phase space is divided in equally- 
sized, five-dimensional cells. For each of these, the expected number of events, ngt, is 
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obtained by the MC simulation described above, after scaling the MC sample to the number 
of events in the data sample. For each cell with a minimum expected occupancy of 5 
events, a is computed comparing the predicted population to the observed number of 
events, (^obs — ""-fit)/*^) where a = y^riohs- The total is formed by summing over all cells 
with minimum occupancy. The number of degrees of freedom is the number of cells with 
minimum occupancy minus the number of free parameters in the fit. The baseline model 
has 9 amplitudes, with 16 free parameters. There are 166 bins with at least the minimum 
occupancy, so the number of degrees of freedom is 150. These bins include 97% of the data. 
The total is 348.0, which corresponds to a confidence level (CL) of 10~^^. We conclude 
from the test that the baseline model does not provide an accurate description of the 
data in the five-dimensional phase space. 

We have looked at the distribution of the throughout the phase space. That might 
reveal the need of a missing resonant amplitude, or problems with the representation of the 
existing amplitudes. We found that there are 33 bins with > 6, containing approximately 
15% of the data events. The bins that are responsible for the poor confidence level of the fit 
are evenly distributed throughout the phase space, instead of being concentrated in specific 
regions. In general, these 33 bins arc surrounded by bins having a good x^- 

Addition of further amplitudes of the type D RR do not improve the fit quality. The 
inclusion of an uniform nonresonant term gives rise to large interference with the (77r+7r~ 
mode without improving the likelihood. A fit using an alternative model, in which the 
D — > (77r"'"7r~" mode is replaced by a uniform nonresonant amplitude, is worse than our 
best fit by more than 200 units of w = — 21og£. No improvement is obtained when other 
amplitudes of the type D Rtt^tt^ are added. The contribution from p^7r~^7r~ is 

negligible. This is also the case for amplitudes with higher mass states, like the p(1450)°7r"'"7r~ 
and the /o(1370)7r+7r-. 

It is very unlikely that the poor CL is caused by problems with the representation of the 
signal amplitudes. The dominant resonance is the p(770)°, which has a very well determined 
fine shape. As we will discuss on Section VIE, the uncertainty in the fine shape of the S- 
wave components does not have a large impact on the fit result. The other ingredients of 
the signal amplitudes arc the standard angular distributions and the well known, widely 
used Blatt-Weisskopf form factors, depending on one single parameter whose value does not 
significantly affect the fit result. 
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The isobar model is a very standard tool, and has been successfully used for several 
decades. Almost all Dalitz plot analyses of heavy flavor decays were performed by rep- 
resenting the signal distribution as a coherent sum of resonant amplitudes with constant 
coefficients. The isobar model has also been successful in describing four-body decays in the 
past 
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3J] (although no CL was quoted in these analyses) and recently [35|, l36|, but in all 



cases the data samples have limited statistics. This analysis is the first attempt to use the 
isobar model for an amplitude analysis of a four-body decay with moderately high statistics. 
The underlying picture of the isobar model may be too simplistic, signaling other types of 
effects should be taken into account. 



B. General features of the best solution 

The dominant contribution to the 7r^7r+7r^7r+ decay comes from the — > 

ai(1260)"'"7r~ mode, accounting for 60% of the decay rate. The second most dominant 
contribution is the p(770)°p(770)° mode, with a decay fraction of 25%. Modes of the 

type D — > Rn^n^ correspond to 20% of the decay rate. 

The large value of the —>■ ai(1260)+7r^ decay fraction supports the picture of an ex- 
ternal VT-radiation amplitude, with the virtual W coupling to the ai(1260), as the dominant 
mechanism for the 7r^7r+7r^7r+ final state. The same dominance has also been observed in 



other four-body decays, such as K tt+tt tt^ and —>■ K^n^n tt"*" 33|, l3J]. The 



same picture can be drawn in the case of three-body decays with the W coupling to the 



p(770)"'', such as —>■ K^-k^tt^ 37|] or —>■ K^ti^'k^ [38], where the contribution from 
D Kp exceeds 50%. This pattern is similar to the well known vector-dominance in elec- 
tromagnetic interactions. It can be understood as a manifestation of the V-A nature of the 
weak interaction. 

The relatively large contribution (25%) of the D p^p^ {D —>■ V1V2 ) decay has also 
been seen in other final states, such as —>■ K'tt'^tc^-k^ , where the mode —>■ K*^p^ 



accounts for 17% of the decay rate, accor ding to [33|| or over 40%, according to [34]. In the 
case of the final state K~ K^tt^tt^ 36|], the contribution of the modes p^cj) and 



K*^K*^ amounts to more than 30%. In all cases the dominant D —>■ V1V2 amplitude 
is the internal W^-radiation, which is expected to be suppressed with respect to the external 
VT-radiation amplitude. 
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The remainder of the decay rate is due to three-body nonresonant modes D f^vr+vr". 
Like the D p^p^ decay, the internal ly-radiation should be the dominant amplitude. For 
any possible meson R there are several ways to form a Rtt^ti^ state with J = 0, combining 
the orbital angular momentum between the two pions from the nonresonant tt+tt^ system, 
the orbital angular momentum between this system and the resonance R, and the spin 
of the resonance. Consequently, for any possible resonance R we can assign different spin 
amplitudes. We tried all possible assignments for the p(770)°7r+7r~ and for the p(1450)'^/T+7r~ 
spin amplitudes, but none yielded a significant contribution. It is interesting to note that 
the only D — > Rti^ti^ modes with a significant contribution are an^Ti^, /o(980)7r+7r^ and 
/2(1270)7r+7r^, with the nonresonant tt+tt^ pair being in a pseudoscalar state. The tt+tt" 
pair and the resonance must then be in a state of even relative orbital angular momentum, 
and the resulting spin amplitudes coincide with the ones from three-body D decays 



C. ai(1260) results 

As we mentioned before, the study of the ai(1260) is one of the main purposes of this 
analysis. The ai line shape is connected to the resonant substructure through Eq. iQ. The 
data seem to require, in addition to the dominant af — > p(770)^tt^ amplitude, an S-wave 



component af —>■ (J7r+, as in the case of CLEO 3l| and E852 39j. A fit without the an^ 
mode increases the quantity w = — 21ogi2 by more than 200. We tried other modes, such 
as the af /o(980)7r+, af p(1450)°7r"*" and af — ^ /o(1370)7r'*', but their contribution 
is negligible. We have also made a fit with the —>■ /o(980)7r+ and a'l /2(1270)7r+ 
amplitudes replacing the /o(980)7r+7r^ and /2(1270)7r+7r^. This fit is worse by more than 
100 units of w = — 21og£. 

We found that the three channels for the ai resonant substructure are necessary and 
sufficient for a description of the data. Our model for the ai substructure is the same as 



3l|. The 



that of E852, but it is in contrast with the CLEO analysis of z/^7r+7r'^7r° 
CLEO fit to the tt^tt^tt^ line shape required seven different contributions. However, the 
three analyses agree on the general picture of the ai substructure, according to which the 
dominant contribution comes from the S-wave p(770)°7r, (over 60% of the ai decay rate), 
followed by the an (~15%) and with a small D-wave p(770)°7r component. 

From the results of the ai resonant substructure we can measure the ratio between the 
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D- and S-wave p°7r+ amplitudes. In the flux-tube-breaking model [40|] this ratio is given by 



where D and S are the D- and S-wave amplitudes. Using this deflnition we measure 

A{ai (p7r)g) ^ _Q ^ q qqq ^ q ^24) 
A{ai ^ (p7r)5) 

where the statistical uncertainty is obtained from the full covariance matrix. 

The ai(1260) line shape was determined by the iterative procedure described above. We 
scanned over the values Sq^ and Fq^ in order to flnd the combination of values that optimizes 
the flt result. The resonant substructure remained stable during this scanning procedure; 
that is, the values of q varied less than 10%. We found that there is some correlation 
between Fq^ and the coefficients Cj. Fits of nearly equivalent quality can be obtained 
using combinations of these parameters along the diagonal of a rectangle, in the \/~^ x Fq^ 
plane, with sides 1230 < < 1270 MeV/c^ and 520 < F^^ < 680 MeV/c^. The optimum 

values are = 1240 MeV/c^ and F^^ = 560 MeV/c^. 

The total width FJJ'qj(s) (Eq. 8) is shown in Fig. [HI as a function of the 7t~^7t~7t~^ mass 
squared. The solid line is the total width, whereas the dashed line represents the contribution 
from the K*K channel. The opening of this channel introduces a cusp-like effect on the total 
width. 

Our values of the ai parameters differ significantly from CLEO result: a/s^ = (1331 ± 
10) MeV/c^ and Fq^ = (814 ±38) MeV/c^. The CLEO values, however, depend strongly on 
their form factor parameter r. When the value of this parameter is set to r = 1.2 GeV~^ in 
the CLEO analysis, their values of \/s^ and Fq^ are in very good agreement with the result 



'^'-'-'^ ^ 

of this analysis. Our result is in good agreement with the model of Kiihn and Santamaria 



4l| {^/^ = (1262 ± 11) MeV/c2 and F;;^ = (621 ± 66) MeV/c^) and a bit higher than the 



values from Isgur et al. = (1210 ± 7) MeV/c^ and Tq' = (457 ± 23) MeV/c^) |40|. 

Both models consider only the pvr mode for the ai resonant substructure. 

D. 1)0 p^p'^ results 

The amplitude for the D — > p^p^ decay is written in the helicity formalism. The most gen- 
eral amplitude D Vi(Ai)V2(A2) is a sum of three independent components, corresponding 
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FIG. 7: The difference -21og>C + 21ogi;max as a function of the Tq^ parameter, for ^/s^ = 1240 
MeV/c2. 
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FIG. 8: The s-dependent oi width, r"^^(s), plotted as a function of the tt+tt tt"*" mass squared, s. 
The dashed line represents the contribution from ai K*K. 
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to the linear polarization states of the vector mesons. The is an initial state with J = 0, 
so the vector mesons helicities Ai and A2 must be equal, Ai = A2 = A = —1, 0, +1. The three 
independent helicity amplitudes, corresponding to each of the possible helicity states are 
H_, Hq, and Alternatively, one can decompose the D ^ V1V2 decay amplitude in the 
so-called transversity basis [42], with amplitudes given by Aq = Hq, A\\ = {H^ + H_)/\^ 
and A_i_ = {H^ — H^)/\/2. This basis represents states in which the vector mesons are 
longitudinally polarized {Aq), or transversely polarized (Ay, A_[_). The transversity basis 
states are parity eigenstates. 

In the D rest frame the two vector mesons are back-to-back. After a boost along the 
line-of-flight to the Vi rest frame, we define the helicity angle 9i as the angle between the vr" 
and the direction of the boost, measured counterclockwise. The helicity angle 62 is defined 
in a similar way. With these definitions and the acoplanarity angle, Xi defined as the angle 
between the Vi, V2 decay planes measured in the D frame, one can write the D V1V2 
decay amplitude as 



A{D — > pp) oc Aq cos 6*1 cos 02 H — J= sin 61 sin 6*2 cos x ^ sin ^1 sin 6*2 sin x- (25) 

V 2 V 2 

In the above equation the amplitudes Aq, A\\, and A^ are of the form Ai = aiC^^^ f {m^^) , 
where Oj and 5i are real constants and f{rn^^) has all the energy dependent terms: Breit- 
Wigner functions, form factors, and the energy dependent part of the angular distribution. 

The longitudinal polarization is defined as the ratio of the longitudinal to the total decay 
rate, 



\A 



(26) 



Using the values obtained by the fit (Table HTB we measured Pl = (71 ± 4 ± zj/o. 

Our value is in agreement with the BaBar result from B'^ — > [l^, Pl = (87 ± 
13 ± 4)%, but is somewhat different from the polarization measured in —>■ p^p^ {Pl = 
(99±3±4)%) |lo| and in B^ ^ p^p^ {Pl = (95±11±2)% [9.] and Pl = (90.5±4.2±2.5)%) 



ll|. 



The polarization measured in B decays is in agreement with theoretical expectations 



which are based on the assumption that the transition amplitude for the B pp can be 
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written as a product of two independent hadronic currents. The factorization hypothesis 
may be a good approximation for B meson decays, due to the large value of its mass, but 
fails to describe hadronic decays of D mesons. That could explain the observed difference 
in the values of the longitudinal polarization fraction from D and B decays. 



E. Systematic uncertainties 

The representation of the structures present in 7i~7t~^7i^7t^ is a major source of systematic 
uncertainty. The problem starts with the choice of a set of amplitudes to describe the data. 
Furthermore, our model depends on some input parameters, which are fixed during the fit: 
the parameter for the form factors; the ai parameters Sq^, Fq^ and gK*K', the relative 
amount of the backgrounds. In addition, there are the parameters defining the a and the 
/o(980) line shape: rriQ and Fq, and Qt^, Qk, and mg". We should also mention that there is 
some freedom in the choice of the angular distributions. 

We investigated the effect of the uncertainty on the input parameters of our model on 
the fractions and phases for the various intermediate channels. For this purpose, the data 
was fit with variations of our baseline model. We varied the value of from 1 to 5 GeV~^. 
We have also fit the data using different sets of ai parameters, s^^ , Fn^ and Qk'K- For the 
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45| 



cr and /o(980) line shapes we used the BES parameterization 

In addition, we varied the background parameterization: the relative amount of back- 
ground, the background composition, adding a possible contribution from uncorrelated 
events, and the relative fractions of the backgrounds, changing the fraction of p^n^Ti^ from 
to 50%. 

A very small change of the values of the fit parameters, compared to their central val- 
ues, is observed when we vary the parameter, the background parameterization and the 
/o(980) line shape. The systematic errors are largely dominated by the uncertainty in the ai 
parameters and on the a line shape, in this order. The total systematic errors are taken as 
the root mean square of fit variations. These are propagated to fractions, the longitudinal 
polarization measurement, and the D/S ratio. The systematic errors are quoted in Table 

Iml 
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VII. CONCLUSIONS 



In this paper we presented a new measurement of the relative branching ratio r(D° — >■ 
7r+7r-7r+7r-)/r(D° K-^+^-tt+) = 0.0914±0.0018±0.0022, based on a veryclean sample 
of 6360±115 events. Our value is compatible with the recent CLEO-c result [20;], which has 
a sample of comparable size, but with a higher background level. 

A full, coherent amplitude analysis of the ti^ti^ti^ti^ decay was made for the first 

time. The isobar model with nine resonant amplitudes provides a description of the general 
features of the data. In our model there are contributions of three types. The dominant 
contribution comes from the decay D° ai(1260)"''7r~, accounting for 60% of the total decay 
rate. This is followed the decay —>■ p(770)°p(770)°, whose relative fraction amounts to 
25%, and by the three-body nonresonant decays an^n^, /o(980)7r+7r^ and 

D° /2(1270)7r+7r-, with a combined fraction of 20%. 

The —>■ TT~7i~^7T~7r~^ decay is a very suitable tool for investigation of the ai(1260) meson. 
A good description of this state is achieved assuming a simple resonant substructure, with 
only three contributions: a dominant p(770)°7r+ S-wave, a small —>■ p{770)^7T~^ 

D-wave, and the cxtt^. In the framework of the flux-tube-breaking model, we measure 



the at p(770)°7r+ D/S ratio to be = -0.043 ± 0.009 ± 0.005. 

The line shape of the ai(1260) is also an output of the amplitude analysis. The best 
description of our data is achieved with the ai(1260) Breit-Wigner parameters, a/s^ and 
Tq\ in the range 1230 < < 1270 MeV/c^ and 520 < Tq' < 680 MeV/c^ with optimum 

values at = 1240 MeV/c^ and T^^^ = 560 MeV/c^. 



The — s> 7r^7r+7r^7r+ decay has an important p^p^ component. The helicity formalism 
is used to to describe the angular distribution of the four pions from this mode. Using 
the transversity basis states and a full amplitude analysis, we measure the ratio of the 
longitudinal polarization to the total D —>■ p^p^ rate to be Pl = (71 ± 4 ± 2)%. 

The fit quality, from an inspection of the five-dimensional phase space, is poor and it is 
not improved by changing the representation of individual signal amplitudes or by adding 
other amplitudes to the baseline model. We conclude that in the case of —>■ 7r^7r+7r^7r+ 
the failure of the isobar model in providing an accurate description of our data suggests that 
other type of effects should be taken into account. Possible candidates are the Bose-Einstein 
correlations or a non-uniform nonresonant amplitude. The most crucial missing ingredient 
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may be the energy dependent final state interactions. Unfortunately this is a very difficult 
problem, even for three-body decays. To our knowledge, there are no calculations of the 
effects of FSI in four-body decays. This problem, however, must be addressed in order to 
explore the full potential of the hadronic decays of heavy flavor and the amplitude analysis 
technique. 
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APPENDIX A: ANGULAR DISTRIBUTIONS 



The angular distributions of the final state particles result from angular momentum 
conservation throughout the whole decay chain. These distributions depend on both the 
spin of the intermediate resonances and on the orbital angular momentum at each vertex. 
In the case of sequential decays, like the D a^vr" , p°7r"^, p° ti^ti^ chain, parity 

conservation is also a constraint on the Oi decay and all its subsequent decay products. 

Amplitudes describing angular distributions, hereafter called spin amplitudes for short, 
are, in general defined using the relativistic tensor formalism 



251, 



26j. This formalism ex- 



plores the connection between the only available observables — the momenta of the final 
state particles — and the spin/orbital angular momentum dynamics. Polarization vectors, 
representing the spin/orbital angular momentum involved at each vertex, as well as all rele- 
vant 4-momenta, are written in terms of these final state momenta. Covariance, invariance 
under rotations and parity conservation, whenever applicable, restricts the many possible 
combinations between 4-momenta and polarization vectors. 

The relativistic tensor formalism was used in this analysis to describe all intermediate 
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channels, except for the P^P^i for which we used the hehcity formahsm. 

1. — > aiTT 

In our model for the ai resonant substructure there are three amplitudes: ai p'^tt^, 
with the p° and the vr"*" in relative S- and D-wave; and the ai air^. We denote the ai 
and p° polarization 4-vectors with helicity A by e^{\) and e^(A), respectively. 

We form Lorentz scalars at each step of the decay sequence by contracting these polar- 
ization vectors with the appropriate 4-momenta. In a step where the resonance appears as a 
decay product, the complex conjugate of its polarization 4- vector is used. If the ai is formed 
by particles 1, 2 and 3, then the amplitude for the a'^n^ is 

(ai7r4 |A| oc £;(AK. (Al) 
The amplitude for the ai p^n^ (p°^^ iii S-wave) decay is 

(p7r|A|ai) oc£,(A)e*"(A'), (A2) 

and (p°7r+ in D-wave) 

(p7r|A|ai)cxe.(A)g^e:(A'K^, (A3) 

with Q =pp- pp, p^=p-^+ p2 and = p^ + p2 + p^ 
The amplitude for the p° decay is 

(ttttIAIp) oce„(A')g", (A4) 
with q = pi — P2- Finally, the amplitude for the an^ is 

(avr |A| ai) oc e^(A)g^ (A5) 

We then combine the amplitudes for each step and must sum over the oi and p° polariza- 
tions, since none of these states is observed directly. The sum over polarizations is simplified 
by using the completeness relation |26| 
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J2 = <(A)£.(A) = -g^, + p^p^^Ip". (A6) 

A 

After summing over the unobserved polarizations, the spin amphtude for — > a^^n^ i 
af — > p^TT"*" S-wave, — > tt'^tt" is 

M=P.(-9.„ + ^j(j --r)^,- (A7) 
For the af — > p^Tr"*" D-wave the spin amphtude reads 

M=p', ^-9,. + ^) g.p^, + ^) (A8) 

FinaUy, the spin amphtude for the af — > air'^ is 

A< = -^mv + Q'- (A9) 



ai 



2. p^p^ 

The spin amphtudes for the D p^p^ decay are written in the hehcity formahsm. This 
decay is of the type P — > Vi(Ai)y2(A2), where Ai, A2 are the vector meson hehcities. Given 
that the initial state has J — M — the constraint |Ai — A2I < M imphes that Ai = A2 = A. 
The only allowed values for the vector meson hehcity are A = —1, 0, and +1. The helicity- 
basis states |JM, A) are |00,-1), |00,0) and |00,+1). 

For each helicity state there is one independent amplitude. These amplitudes are func- 
tions of the helicity angles, defined in each resonance rest frame, and of the angle formed 
by the vector mesons decay planes (acoplanarity angle), measured in the D rest frame. We 
define the z axis along the line of flight of the vector mesons in the D rest frame. With the 
notation p° T^iT^t ~^ '^^'^t ^ helicity angles 61 and 62 are defined as the angle 

between tti and the z axis, in the p\ rest frame, and the angle between tts and the z axis, in 
the P2 I'sst frame. 

Now with all momenta measured in the D frame, the orientations of the decay planes are 
n\ — p(7rj) X p(7rf) and 712 — p{T^t) x p{t^^)- The acoplanarity angle is x = fii • 77,2. 

With the above definitions of 9i, 62 and x, the general form of the hehcity amphtudes is 
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A, = e^'^dlMdlM. (AlO) 
With the appropriate d functions we have 

A_i = ^e"^^ sin 61 sin 62, (All) 

^0 = cos 6*1 cos 6^2, (A12) 



A+i = -e^^sin^isin^2- (A13) 

The hehcity-basis states are not parity eigenstates. One can construct a new basis with 
definite parity, 

l/„ >^ ^J''^-'\ (A14) 



V2 



|00,1) - |00,-1) 
\h >= ^ , (A15) 



|/l>=|00,0). (A16) 

n 

This is known as the transversity basis [4^]. The corresponding amphtudes are 

A\\ = "^"^ = cos Ysin6'i sin6'2, (A17) 

A± = — — = —= sin X sin 9i sin 62-, (A-18) 

V 2 V 2 



Al = Ao = cos Oi cos 02. (A19) 

The above amphtudes are expressed in terms of non-covariant quantities, defined in dif- 
ferent reference frames. Covariance is recovered by muhiplying the hehcity amphtudes by 
an energy dependent term, which accounts for the boosts from each resonance frame to the 
D frame. Each of the above amphtudes is then multiphed by 
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^a-^axJ^^, (A20) 
where s is the two-pion invariant mass, p*{s) is the break-up momentum and r(s) is 

These are the formulae used for the D p^p^ spin amphtudes. 
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